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I OBJECTIVE (U) 


, The objective of this effort is to investigate the 

possible effects of ambient geophysical/ low-frequency electromagnetic 

' ‘ > . 

factors on remote viewing (RV)* performance 


4 


(U) RV (remote viewing) is the acquisition and description, by mental 
means, of information blocked from ordinary perception by distance 
or shielding. 
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II INTRODUCTION (U) 

K SRI International is tasked ^ 

'to investigate a potential correlation between remote viewing 
(RV) performance and ambient geophysical/extremely-low-frequency electro- 
magnetic (ELF) activity. The possibility of such correlation is indicated, 
for example, by studies showing psychophysiological effects 1 ' 3 * and 
behavioral changes ' associated with ELF electro-magnetic fields. The 
geophysical variables of interest include such factors as ELF intensity/ 
fluctuations, ionospheric conditions, geomagnetic indices, sunspot number, 
and solar 'emissions (e.g., X rays and solar flares). The questions to be 
answered in. this program are 

• Do geophysical/performance correlations exist such that 
measurement of the ambient geophysical variables could 
be used as an indicator of expected performance? 

• If so, can optimum performance windows be identified? 

(U) The structure of the program that will address the above issues 
consists of 

• A literature search. 

• Real-time ELF measurements 

— SRI (Menlo Park, California location) 

— Time Research Institute (Los Altos, California 
field station). 

• Real-time geophysical data acquisition via NOAA 
(National Oceanic and Atmospheric Administration) 

Westar IV satellite downlink. 

• Computer correlation studies of RV performance versus 
variables of interest. 


(U) References are listed at the end of this report. 
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III PRESENT STATUS (U) 


A. (U) General 

(U) In order to accomplish the goals set out in Section II, the 
program has been designed as a joint effort between SRI International and 
Time Research Institute of Los Altos, California, with SRI as the prime 
contractor. Time Research Institute is a research organization that 
specializes in temporal analysis of geophysical variables and their po- 
tential correlation with phenomena of interest, such as weather patterns, 
earthquakes, and so forth. 

(U) With regard to the present effort, Time Research Institute is 
responsible for establishing the appropriate hardware and software systems 
for collecting and analyzing data concerning environmental conditions and 
their correlation with RV performance. The purpose of the correlation 
study is to determine whether RV performance is enhanced or degraded by 
measurable changes occurring in the geophysical (including solar- 
terrestrial) environments. The specific data bases under consideration 
in this effort are given in Table 1. 


(U) Should correlations between geophysical variables and RV 
performance be found, the application potential of the effort is twofold: 

(1) Time periods in which enhanced RV performance might be 
expected could be identified, resulting in increased 
quality and accuracy of information obtained through 
such channels; similarly, time periods in which degraded 
RV performance might be expected could be avoided. Thus, 
optimum performance windows would be identified. 

(2) An increased understanding of the types of environmental 
changes that correlate with RV performance could provide 
clues as to the mechanisms involved in RV functioning. 

Such knowledge would lead to more focussed research on 
factors that could enhance RV performance, and would 
also provide information critical to the development of 
defensive countermeasures against RV. 
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Table 1 

(U) GEOPHYSICAL DATA BASES 

• Solar-Terrestrial 
— Geomagnetic 

Ground-measured indices Ap, K, etc* 
Satellite-measured intensities 

— Solar flux (MHz) 

15,400 1,415 

8,800 606 

4,995 410 

2,800 245 

2,695 

— Sunspot number 
— Solar flares 

— Interplanetary magnetic field 
— Solar wind (Pioneer XII) 

— Protons 

“ Cosmic ray indices (neutron monitor) 

• Ionospheric Measurements I 

— Sudden ionospheric disturbances (SIDS) ! 

— Auroral electrojet 

— Radio propagation quality indices 

• ULF/ELF 

— 30 frequencies (from 1 to 30 Hz) 

UNCLASSIFIED 


(U) The tasks listed in Section II (literature search, real-time 
ELF measurements, real-time geophysical data acquisition via satellite 
downlink, and correlation studies of RV performance versus geophysical 
variables of interest) have been prioritized with the goal of producing 
the longest possible ELF data base during the period of this contract* 
Therefore, while all of the tasks are being pursued in parallel, ELF- 
related tasks have been the focus of attention to date. 

6 
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B. (U) Literature Search 

(U) The purpose of the literature search is to acquire and 

integrate information from recent work done in the field of biological 
response to ELF and geophysical conditions. Approximately 30 keywords 
are in use as input to a computerized literature search. Forty manu- 
scripts have been obtained to date and are in the process of being 
reviewed, plus additional sources of literature have been identified and 
will be retrieved as priorities permit. 

C. (U) ELF Measurements 

1* (U) Introduction 

(U) Although the ELF frequency range (3 to 300 Hz) has been 
studied in some detail, many unknowns remain. For example, although it 
is known that ELF frequencies generated by geophysical means (e.g., elec- 
trical storm activity) tend to distribute themselves globally, little 
information is available on the variation of the ELF environment from 
location to location. Therefore, local variations may exist that are 
caused by both manmade sources, and by the geological structure of the 
area. In the San Francisco Bay Area, manmade sources that generate ELF 
on a local scale include motors, telephone lines, power lines, and 
electrical subways [Bay Area Rapid Transit (BART)], and it needs to be 
determined whether the emission from such sources constitutes a signif- 
icant contribution to the omnipresent global ELF field. 

(U) In order to address the above issue, two ELF monitoring 
stations are being set up— one at SRI Menlo Park (in the RV Laboratory), 
the other at the Time Research Institute field station, 17 km distant. 

It is anticipated that the SRI environment may be a "noisy" one due to 
the large amount of electronics known to be in the area. Data from the 
two sites, taking the field site as a reference, are to be compared in 
order to begin to differentiate the naturally-occurring ELF from the 
manmade noise occurring at the location where RV is being carried out. 
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2. (U) Los Altos Site 

(U) Since May 1982, Time Research Institute hqs been operating 
an ELF monitoring site in Los Altos, collecting data twice daily for the 
purpose of correlating ELF disturbances with various phenomena of interest. 
In this period, analysis techniques were developed that are directly 
applicable to the present task. 

(U) One of the first tasks was the upgrading of the Los Altos 
ELF monitoring site to provide coverage during power interrupts. Research 
was performed to determine the best power- interrupt system, including 
generators. The selection criteria chosen for the generator sought to 
optimize power output, cost effectiveness, and reputation for reliability, 
with the consideration that this system might serve as a model for addi- 
tional sites in the client community. On this basis, a 3500-W Kubota 
generator was bought and installed at the site. A PTI Industries "Data- 
shield" device was also purchased and installed, for use in conjunction 
with the generator (which must be started manually) • This device powers 
the required electronics for a period of twenty minutes on its own while 
awaiting generator startup. Furthermore, an automatic alarm telephone 
dialing system dials up as many as four individuals should a power failure 
occur while no one is at the site. The two devices working together have 
protected the system on numerous occasions from power interruptions-- 
including an 8-hr outage planned by Pacific Gas and Electric Company. 

Thus, since installation of the power- interrupt equipment, there has been 
no loss of data collection/storage. 

3 . (U) ELF Data Acquisition Systems 

a. (U) Basic System Design 

(U) With the requirement that two ELF monitoring sites 
be implemented for the program--one at SRI and one at Time Research 
Institute — it was decided that the two systems would be made identical. 

In this way, differences between the two systems would be minimized, thus 
reducing the opportunity for artifactual differences between the two sys- 
tem outputs. 
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(U) Figure 1 is a block diagram of the basic ELF data 
acquisition system. The ELF signal is collected by an antenna, amplified, 
and then digitized by an ana log- to-digita 1 (A/D) converter so that the 
signal can be input into an IBM PC Jr. computer for the purpose of analy- 
sis by a fast-f ourier- transform (FFT) program on at least an hourly basis. 
The recorded data are then transferred by floppy diskette to an IBM XT 
computer for further handling. 



INCREASES SIGNAL STRENGTH 


MAKES SIGNAL READABLE BY 
COMPUTER 


DETERMINES INTENSITY OF SIGNAL 
FREQUENCY COMPONENTS 


FLOPPY DISKETTES TRANSPORTED 
TO IBM XT COMPUTER 


UNCLASSIFIED 


FIGURE 1 (U) ELF DATA-ACQUISITION SYSTEM 


(U) As indicated in the above system description, an 
integral part of the data acquisition system is computerized record 
keeping, using IBM sy stems--bo th the IBM PC Jr. and the IBM XT. The 
software is written and compiled on the XT, used as a master system, and 
then run on the PC Jr. (the PC Jr. is not itself large enough to run a 
compiler, nor are there compilers written for it). Beyond this, however, 
in spite of the much publicized n compa tibility" between the various 
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IBM PC systems, some development time had to be expended to ensure that 
programs compiled on the XT (that were concerned with communication with 
external devices) could be run on the PC Jr. --because they handle the 
addressing of communication ports differently. 

b. (U) FFT Program 

(U) A consultant was hired from the Radio science Labora- 
tory at Stanford University to critique the ELF system at Time Research 
Institute as it stood at program start. As part of that critique, he 
recommended changes in the FFT software to increase its accuracy and 
running speed. As a result, new FFT software was written in compiled 
BASIC and then debugged. The program description is given in Appendix A. 
This task was completed in mid-April. 

c. (U) System Electronics 

(U) The prototype system of Figure 1 was assembled and 
installed for our purposes at the Los Altos site. Early tests indicated 
that modification of the existing electronics was required because 
(1) the system was sensitive to ground-loop problems, thus the preamplifier 
had to be redesigned to include an isolation amplifier, and (2) when the 
new system was installed in mid-June, it was found that a slowly-varying 
dc level was superimposed on the incoming ELF signal, resulting in ex- 
cursions that exceeded the limits of the A/D converter at amplification 
levels required for good signal analysis. Therefore, new circuitry was 
designed to eliminate the dc problem. With these changes, the ELF detec- 
tion system is scheduled for completion of testing, debugging, and 
calibration at the Los Altos site in August, before its sister system 
is installed at SRI. 

d* (U) ELF Electronics/Sof tware Subtasks (Status) 

(U) The status of the ELF electronics/software subtasks 

is as follows: 
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• Subject to the requirement that the basic microprocessor 
units to be used in the program would (1) be able to com- 
municate with an IBM XT, (2) have at least 64K of memory, 

(3) possess diskette storage capability, and (4) be cost 
effective, the IBM PC Jr. was selected from among the 
various alternatives, and two units were purchased. 

• Preamplifiers, low bandpass filters, and amplifiers have 
been designed, and one complete prototype system has been 
assembled. Amplifiers and filters have been assembled 
for the second system to be installed at SRI. 

• Design specifications have been completed for the ELF 
antenna. Assembly of this antenna is pending the results 
of a calibration task (described in Section III.D.3.e). 

• A survey of A/D converters that would be compatible with 
the other system components was completed, and the 
selected units were purchased. 

• Communication between the IBM PC Jr. and the A/D con- 
verter has been established, enabling the computer to 
read the incoming signal. 

• Software has been developed and debugged that: (1) reads 

the communications RS 232 port of the PC Jr. input from 
the A/D converter, (2) performs FFT analysis of the signal, 
and then (3) writes half-hourly averages of 19 different 
frequencies (from 1 to 29 Hz) to a computer diskette. 

Further software refinements will continue to be made, such 
as determining daily maximum values for each frequency. 

This software has been implemented and debugged. The sys- 
tem is ready to begin data acquisition upon implementation 
of the dc- level-controlling hardware. Minor enhancements 
of the software will continue. 

• The Los Altos system is now in operation and is being 
tested. It is anticipated that the SRI system will be 
implemented in August. 


e - (U) ELF System Calibration 

(U) System calibration is proceeding. The Stanford 
University consultant mentioned earlier is an expert in the areas of ELF 
and VLF measurement, antenna design, and spectral analysis. His cali- 
bration program is being carried out at the present Los Altos site using 
specialized, sophisticated instruments from the Stanford Radioscience 
Laboratory. 
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(U) As part of the calibration process, certain power- 
line and power-supply "noise" sources were identified that could have 
produced artifacts in the historical and ongoing data. Therefore, in- 
creased electronic-filtering measures are being incorporated into the 
system, including software changes to recalibrate past data already in 
storage. 

(U) Finally, as part of the calibration process, three 
systems are to be monitored simultaneously: (1) the original unmodified 

Los Altos system on which the historical ELF data have been acquired; 

(2) a new system (described in this report) using an ELF wire-coil 
antenna designed by and on loan from the Radioscience Laboratory at 
Stanford University; (3) the new system, but with a bioantenna (oak tree) 
as an antenna, which is a procedure recommended by Stanford (see Nature 
reprint, Appendix B) . The results from the calibration tests will be 
used to calibrate the historical data base, and to fix the final antenna 
design. 

U« (U) Satellite Downlink Geophysical Data-Acquisition System 

| " 

(U) A near real-time satellite downlink system for $olar- terrestrial 
data has recently become available from NOAA (National Oceanic and Atmo- 
spheric Administration). With this unit, it is possible to provide 
immediate feedback and/or analysis in conjunction with RV: sessions. 
(Normally, there are long delays in procuring solar- terrestrial data; 
without the downlink, delays of 10 days to 6 months are standard.) The 
downlink system provides for accumulation of a detailed data base directly 
on computer diskettes. (See Appendix C for an item-by-item description.) 

(U) A satellite controller and a dish antenna for the downlink 
system were ordered and installed at the Los Altos site early in the 
project (see Figure 2). At the time of this writing, specifications for 
data-acquisition software for the IBM PC Jr. have been completed, and 
first-order software has been written that captures the data to computer 
diskettes. Because of the large volume of data transmitted each minute, 
a double-density, double-sided diskette fills in about 21 hours. Disk- 
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ANTENNA 




UNCLASSIFIED 

FIGURE 2 (U) REAL-TIME GEOPHYSICAL DATA 

ACQUISITION VIA WESTAR IV 
DOWNLINK 


(U) 

ettes have been changed on an almost daily basis since early April in 
order to begin to acquire a long-term data base. 

(U) To increase the amount of data that can be stored on a 
diskette, it is necessary to separate numerical data of interest. The 
necessary software is now being written that will be able to identify 
the various data types as they are transmitted, so that only the data of 
interest will be transferred to an appropriate data file on the diskette. 
Statistical analysis can then be done on the data in this form. 

(U) Other data bases are continuing to be maintained for this 
project by Time Research Institute. Files of 2800-MHz solar flux, the 
planetary magnetic activity index (Ap) , the Anchorage magnetic index, 

13 
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(u) 

and the Stanford mean solar magnetic field are updated on a weekly or 
monthly basis as the data are available, ELF data from the "old" system 
are recorded twice daily; ELF data from the "new 11 system are being 
recorded at the following intervals: half-hourly averages (48 times a day), 
and two sets of half-daily averages (twice a day) at 00:00 and 12:00 UT, 
and at 04:30 and 16:30 UT. Whole-day averages are also being recorded. 

(U) To summarize the status of the satellite-downlink data- 
acquisition system: 

• The downlink-geophysical system is in place and in operation 
at the Los Altos site. 

• Partial data acquisition from the downlink is in place. 

• Software for the final data-acquisition system will be 
completed and implemented in August. 

E. (U) Geophysical Data/RV Performance Correlation Analysis 

(U) Data for RV sessions are to be analyzed for statistical cor- 
relation with respect to the battery of geophysical data sets listed in 
Table 1, and those acquired by Time Research Institute via measurement 
(Figure 3). Recording of RV and geophysical data is now in progress. 

The overall system for data acquisition and analysis of RV performance/ 
geophysical data is depicted in Figure 4. When enough data have been 
collected toward the end of the contract period, analysis will be performed. 
The tasks described in earlier sections are in preparation for this task, 
and therefore have received the bulk of the effort. Certain subtasks in 
the analysis task, which require a longer leadtime, have^ however, already 
been completed in preparation for the analysis. 

(U) The primary statistical program that will be used to scan the 
data for possible relationships is called EPOCH ANALYSIS. This program 
reads two files simultaneously. The first file is an event file, the 
second a data file. The program first reads an event, then scans the 
data temporally backward and forward in time around the event. This 
information is stored, a second event is read, and so forth. When all 
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LOS ALTOS STATION 


SRI STATION 



UNCLASSIFIED 


FIGURE 3 (U) REAL-TIME GEOPHYSICAL DATA-ACQUISITION SYSTEM 


(U) 

the events and surrounding data have been read, a printout is created 
that lists appropriate cross-correlation statistics between event and 
data elements. 

(U) Preliminary scans of data generated during an approximate 
100-site series with one remote viewer have been carried out. The 
session quality was graded on a scale of 0 to 3 + , and correlations 
between solar magnetic field and solar sunspot number were investigated. 
Some correlation between RV performance and solar sunspot number was 
found, which, if substantiated by further data, would indicate the 
possibility that performance might improve immediately after a peak in 
the sunspot number, and would deteriorate just before the sunspot number 
peaks in its 27-to-29-day cycle. This result is based on data points 
that are too small in number to be taken seriously at the point, however, 
and is mentioned only to give an example of the types of correlations 
that will be sought out and examined during the course of study. 
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FIGURE 4 (U) GEOPHYSICAL/PERFORMANCE DATA-ANALYSIS SYSTEM 
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IV SUMMARY (U) 


(U) Approximately 80 percent of the 
tasks have been completed, and all of the 
been purchased and delivered. 


project's da ta-acquisition 
equipment and hardware have 


(U) System calibration should be completed in the near future, 

and ELF and downlink-data acquisition will have begun in their final 
formats . 


(U) At the above point, the focus of effort will turn to analysis 
of past and present geophysical and ELF data, soon to be followed by the 
initiation of correlation studies of these data against RV performance. 
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Appendix A 

FAST FOURIER TRANSFORM ROUTINE FOR ELF DATA (U) 
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Appendix A 


FAST FOURIER ROUTINE FOR ELF DATA (U) 


* 

PROGRAM DESCRIPTION * 


(THIS PROGRAM USES THE FAST FOURIER TRANSFORM ( FFT ) ALGORITHM TO 
{CALCULATE THE SPECTRAL MAGNITUDES OF AN ARRAY OF CLOSE-PACKED REAL INPUT 
{DATA POINTS. THE PROGRAM IS ORGANIZED AS FOLLOWS: 

( 1. AT THE START, DATA POINTS ARE STORED IN THE REAL ARRAY X(1), 

; WHERE 1 RUNS FROM Z TO NX-1. NX MUST BE AN INTEGRAL POWER OF 2, 

; OF THE FORM NX = 2-^MX (IE, MX ■ LOG2(NX>>. SUCCESSIVE ELEMENTS 

{ OF X REPRESENT SUCCESSIVE SAMPLES OF AN INPUT SIGNAL, SAMPLED AT 

{ REGULAR INTERVALS OF TIME DT. 

{ 2. THE REAL ARRAY X< 1> IS TREATED FOR THE FFT AS A COMPLEX ARRAY 

; OF NX/2 PAIRS OF REAL AND IMAGINARY ELEMENTS. THAT IS, THE REAL 

; ELEMENTS OF THE ARRAY ARE THE EVEN-NUMBERED INPUT SAMPLES, AND 

{ THE IMAGINARY ELEMENTS OF THE ARRAY ARE THE ODD-NUMBERED INPUT 

» SAMPLES. AFTER THE FFT IS CALCULATED, AN ADDITIONAL STEP IS USED 

; TO EXTRACT THE SPECTRUM OF THE REAL INPUT DATA. TRANSFORMING 

5 CLOSE-PACKED DATA IN THIS WAY, EVEN THOUGH IT REQUIRES AN 

; ADDITIONAL STEP FOR THE REAL TRANSFORM EXTRACTION, IS FASTER THAN 

s TRANSFORMING A COMPLEX ARRAY OF NX REAL-IMAGINARY PAIRS (2*NX 

? ELEMENTS), WHERE ALL OF THE INITIAL IMAGINARY VALUES ARE ZERO. 

; 3. THE FFT IS CALCULATED AS FOLLOWS: 

{ A. THE INPUT PAIRS IN X ARE REARRANGED IN BIT-REVERSED 

! ORDER. 

{ B. THE FFT IS CALCULATED USING DECIMATION IN TIME, WITH 

; ANGLE ARGUMENTS IN EACH SUB-DFT APPEARING IN NATURAL 

{ (IE, INCREASING) ORDER. 

; C. AFTER THE LAST PASS THROUGH THE FFT ALGORITHM, THE 

5 ARRAY X CONTAINS SPECTRAL VALUES IN NORMAL ORDER, WITH 

{ EACH EVEN POINT A REAL VALUE AND EACH SUBSEQUENT ODD 

! POINT THE CORRESPONDING IMAGINARY VALUE. 

I 4. FOLLOWING THE FFT THE REAL TRANSFORM IS EXTRACTED. NOTE THAT 

1 IF THE SAMPLING INTERVAL (THE TIME BETWEEN SUCCESSIVE SAMPLES) IS 

{ DT, THEN THE TOTAL SIGNAL INTERVAL PROCESSED IS NX*DT (THE FFT 

{ ASSUMES THAT THE MISSING RIGHT END POINT IS THE SAME AS THE FIRST 

{ POINT). THE FFT GENERATES A SPECTRUM CONTAINING VALUES AT 

{ INCREMENTS OF DF IN FREQUENCY, WHERE DF - 1/(NX*DT). IF WE WERE 

? TO TRANSFORM AN ARRAY OF NX COMPLEX DATA POINTS (WHOSE IMAGINARY 

{ VALUES WERE ZERO, SINCE WE ARE CONCERNED WITH A REAL SIGNAL) WE 

{ WOULD GENERATE NX COMPLEX SPECTRAL POINTS RANGING IN FREQUENCY 

; FROM Z TO ( NX- 1 )*DF Hz. HOWEVER, THE POINTS FROM (NX/2+l)*DF Hz 

5 TO ( NX-1 >*DF Hz ARE MERELY THE COMPLEX CONJUGATES OF THE LOWER 

; POINTS, AND CONTAIN NO ADDITIONAL INFORMATION. THIS IS BECAUSE 

; THE NYQUIST RATE, OR THE HIGHEST UNALIASED FREQUENCY THAT CAN BE 

» SAMPLED, IS 1 /2*NX*DF . WHEN WE TRANSFORM NX REAL CLOSE-PACKED 

{ POINTS AND THEN EXTRACT THE REAL TRANSFORM, WE GET ONLY THE LOWER 

{ NX/2 COMPLEX SPECTRAL VALUES, WHICH ARE ALL THAT ARE NEEDED. 
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5. NEXT, THE SPECTRUM IS CONVOLVED WITH A SHORT WINDOW FUNCTION. 
THE REAL AND IMAGINARY < EVEN AND ODD) POINTS ARE CONVbLVED 
SEPARATELY, SINCE THE WINDOW FUNCTION IS A SEQUENCE OF REAL 
NUMBERS. WINDOWING IS NECESSARY IF THE SIGNAL FILTERS WE ARE 
EFFECTIVELY SYNTHESIZING ARE TO HAVE A USEFUL SHAPE. WITHOUT 
WINDOWING WE WOULD FIND THAT EACH SPECTRAL FILTER WOULD HAVE A 
NARROW PASSBAND BUT SIGNIFICANT SIDELOBE RESPONSES. THAT IS, THE 
MAGNITUDE OF A SPECTRAL LINE X< 1 ) WOULD DEPEND NOT ONLY ON SIGNAL 
COMPONENTS NEAR 1*DF IN FREQUENCY, BUT ALSO SIGNIFICANTLY ON 
COMPONENTS AT OTHER FREQUENCIES AS WELL. WINDOWING BROADENS THE 
SHAPE OF THE PASSBAND OF EACH SPECTRAL FILTER BUT DECREASES THE 
SIDELOBE RESPONSES. THE AMOUNT OF BROADENING AND SUPPRESSION OF 
SIDELOBES DEPENDS ON THE WINDOW ORDER, OR THE LENGTH OF THE 
WINDOW FUNCTION WHICH IS CONVOLVED WITH THE RAW SPECTRUM. 

6. FINALLY, THE MAGNITUDE OF THE SIGNAL AT EACH SPECTRAL 
FREQUENCY IS CALCULATED AS THE ROOT SUM OF THE SQUARES OF THE 
REAL AND IMAGINARY SPECTRAL COMPONENTS (IE, THE VECTOR MAGNITUDE 
OF EACH COMPLEX SPECTRAL POINT) AND SCALED TO MAKE THE PROCESSING 
GAIN INDEPENDENT OF THE SIZE OF THE BLOCK LENGTH NX. WE COULD 
ALSO CALCULATE THE PHASE OF THE SPECTRAL POINTS BUT TiHIS 
INFORMATION ISN'T TOO VALUABLE FOR OUR USE. NOTE THApT THE 
HIGHEST MEANINGFUL SPECTRAL FREQUENCY DEPENDS ON THE ANTI- 
ALIASING LOW-PASS FILTER THAT IS USED WHEN THE SIGNAL! IS SAMPLED. 
THERE IS NOT MUCH MEANING TO SPECTRAL POINTS ABOVE THE FILTER 
CUTOFF FREQUENCY, WHERE SIGNALS ARE ATTENUATED AND FREQUENCY 
ALIASING BECOMES A PROBLEM. THUS, IN THIS PROGRAM WEI DO NOT 
CALCULATE SPECTRAL COMPONENTS ABOVE ABOUT 30 Hz (THE NOMINAL 
CUTOFF FREQUENCY OF THE LOW-PASS FILTER). 

THIS PROGRAM IS WRITTEN IN BASIC FOR THE IBM PC COMPUTER. HOWEVER, IT 
IS EASILY ADAPTED TO OTHER MACHINES. FOUR THINGS TO WATCH OUT FOR WHEN 
TRANSFERRING THE PROGRAM TO ANOTHER VERSION OF BASIC ARE: 

1. ALL ARRAYS START WITH INDEX 0. THAT IS, THE ARRAY X(NX> RUNS 
FROM X ( 0 ) TO X(NX-l). ARRAY REFERENCES WILL NEED TO BE CHANGED 
IF THE PROGRAM IS TO BE USED ON MACHINES WHERE THE FIRST ELEMENT 
OF ARRAYS HAS INDEX 1. 

2. ALL VARIABLES WHOSE NAMES START WITH THE LETTERS "I" THROUGH 
"N" ARE IMPLICITLY INTEGERS (DEFINT I-N STATEMENT). THESE 
VARIABLES MAY NEED TO BE RENAMED. IT IS IMPORTANT THAT LOOP 
COUNTERS AND ARRAY INDEX VARIABLES BE INTEGERS BECAUSE INTEGER 
ARITHMETIC (WHERE AVAILABLE) IS FASTER THAN FLOATI NG-lPOI NT 
ARITHMETIC. 

3. THIS PROGRAM USES THE INTEGER DIVISION OPERATOR "V. THIS 
MAY BE CHANGED TO H /" IN OTHER VERSIONS OF BASIC. 

4. THIS PROGRAM USES LOGICAL OPERATORS ("NOT", "AND", B OR" ) ON 
INTEGER VARIABLES (NOT LOGICAL VARIABLES) IN THE BIT-REVERSAL 
ROUTINE. IF THESE OPERATORS ARE NOT AVAILABLE THE ALTERNATE 
VERSION OF STATEMENTS 3000-3190 USING ONLY INTEGER ARITHMETIC 
MUST BE SUBSTITUTED. 
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UNCLASSIFIED 


; * 
»* 

I** 


CHOICE OF SAMPLING TIME AND BLOCK SIZE 


* 

* 

Hr 


;THE SAMPLING TIME DT AND THE BLOCK SIZE NX TOGETHER HFTFRMTnf tut 
{NUMBER OF SPECTRAL POINTS CALCULATED AND THEIRSPACING IN FREQUENCY 
’THE CHOICE OF SAMPLING TIME ALSO AFFECTS THE FREQUENCIES AT^WHICH^POUFR 
' J, N L2 A . R ,?? NIC INTERp ERENCE THAT IS PASSED BY THE ANTIALIASING LOW-MS?" 
i TisI E Su WILL APPEAR IN THE SPECTRUM. THE FOLLOWING DISCUSSION ASSUMf£ SS 
{THAT THE 30 Hz LOW-PASS FILTER WILL BE USED TO ATTENUATE HIGH fpfonFNrv 
{SIGNAL COMPONENTS BEFORE SAMPLING. A tNUATc HIGH FREQUENCY 

{ 

{IF THE INPUT SIGNAL THAT IS SAMPLED CONTAINED ONLY FREQUENCIES rfi nu 
,30 Hz IT WOULD BE SUFFICIENT TO SAMPLE AT THE NYQUIST RATE JfbS SAMP, F<{ 
{PER SECOND. HOWEVER. SINCE THE FILTER DOES NOT HAVE INFINITE SAMPLES 
’ A JT NATION ABOVE 30 Hz IT IS NECESSARY TO SAMPLE AT A SOMEWHAT HIGHFH 
{30 T Hz. AND THEN DISCARD those SPECTRAL POINTS REPRESENTING SIGNALS ABOVE 
, 

5 I;-rr 0NE CONSIDERATION IN CHOOSING THE SAMPLING RATE IS THE ACTIIAI 

, RATE AVAILABLE WITH A GIVEN A/D CONVERTER. WITH THE CMC BUSSTER^ 

' A ^ CONVERTER, SAMPLING TIMES CAN BE AS SMALL AS 0 0005 s 

HOWEVER, IF MORE THAN ONE SIGNAL IS TO BE DIGITIZED AT ONE TIME 
THE MINIMUM SAMPLING TIME IS 0.01 s PER CHANNEL OR A SAMP I * 
RATE OF 100 SAMPLES/SECOND. THIS SEEMS LIKE ^ReJsoSaB^cJSkE . 

2. A SECOND CONSIDERATION IS THE EFFECT OF THE SAMPI tnc patc nu 
THE FREQUENCIES OF ALIASED POWER LINE HARMONICS ONLY HARMONICS 
AT 60 AND 120 Hz ARE LIKELY TO BE A PROBLEM WITH DT - 0 s 
THE NYQUIST RATE IS 50 Hz, SO THE SPECTRUM WILL CONTAIN PO?ntI’at 
FREQUENCIES FROM 0 TO 50 Hz. ANY 60 Hz SIGNAL THAT IS d! 2 ?I?Fn 
WILL APPEAR IN THE SPECTRUM AT 40 Hz, WHICH POINT WILL BE THROWN 
OUT, SO 60 Hz INTERFERENCE WON'T BE A PROBLEM HOWEVER 

tKII R fI5t N S5sT T b| 2 KE?^ Y; L Ml575SE% A I N JfYZ?NG\;E HE Dl?r iR (!;iGJER 

ara's s:: e js d a Io at 24 *■ «"• ~ « i's™* 

™ e s.?fS'Ksi & ;s e e ^;s?l g ? pac, " s 


NX 

64 

128 

256 

512 

1024 

2048 


DF 

1.5625 Hz 
0.7813 Hz 
0.3906 Hz 
0.1953 Hz 
0.0977 Hz 
0.0488 Hz 


NUMBER OF POINTS <« 30 Hz 

19 

38 

76 

153 

307 

614 
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WINDOWING 


* 

* 

* 


! after the spectrum of the input signal has been found, it is convolved 

iWITH A SHORT WINDOW SEQUENCE TO IMPROVE THE SHAPE OF THE SYNTHESIZED 
FILTERS THOUGH THE SPECTRAL VALUES ARE COMPLEX, THE WINDOW SEQUENCE IS 
A SEQUENCE OF REAL NUMBERS, SO THE CONVOLUTION IS DONE SEPARATELY FOR 
!tHE REAL AND IMAGINARY IeQUENCES IN THE SPECTRUM. THE CONVOLUTION IS OF 
j THE FORM 


t 


X(j> <- W3*IX< j-3)+X( j+3)I +W2*CX< j-2)+X< j+2>3 
+W1 MX(j- 1 >+X( j + 1 >J +W0*X(j) 


!uhfrf y (■,) IS THE REAL OR IMAGINARY SUBSEQUENCE OF THE SPECTRUM, AND THE 
CONVOLUTION IS ILLUSTRATED FOR A 4TH ORDER WINDOW SEQUENCE.; NOTE THAT 
tIF THE CONVOLUTION IS TO BE DONE IN-PLACE, A SMALL ARRAY MU$T BE USED TO 
H^O THE PREVIOUS VALUES OF X(j-3), X(j-2). X(j-I>, and X(j> FOR USE IN 
l SUBSEQUENT CONVOLUTIONS. 

!thf VALUES TO USE FOR DATA POINTS BEYOND THE ENDS OF THE SPECTRUM ARE 
J5| COMPLEX CONJUGATES OF THE POINTS REFLECTED ABOUT THE ENDS, SO THE 
jAUXILLIARY ARRAY OF PREVIOUS VALUES MUST BE INITIALIZED AS FOLLOWS. 

I Re<-3) - Re( 3 > Im(-3> 

, Re( -2 > - Re(2> Im(-2> 

, R®( - 1 ) * Re< 1 ) Im( -l ) 

, PO I NTS AFTER THE FAR END OF THE SPECTRUM ARE REFLECTED IN THE SAME WAY, 
BUT SIN« WE WON ' T CALCULATE SPECTRAL POINTS ABOVE 30 Hz THERE IS NO 
tNEED TO WINDOW POINTS THERE. 

tTHE FOLLOWING WINDOW SEQUENCES ARE THE MOST IMPORTANT FOR OUR USE. 

THESE ARE ALL MINIMUM-SIDELOBE WINDOWS WHICH SUPPRESS SIDELpBE RESPONSES 
ASMUCH AS POSSIBLE FOR A GIVEN WINDOW ORDER, AND HAVE ^EEN; NORMALIZED 
FOR UNITY GAIN. NOTE THAT THE ORDER-1 WINDOW CORRESPONDS TO NO WINDOW 
!a? ALL. THE ORDER-2 WINDOW IS SIMILAR TO THE HAMMING FUNCTION. 


-Im( 3 > 
- Im< 2 ) 
- Im( 1 ) 



ORDER 

I 

2 

3 

4 


W0 

W1 

W2 

W3 

1 .00000 
0 
0 
0 

1 .00000 
-0.42875 
0 
0 

1 .0000000 
-0.5895613 
0.0922278 
0 

1 .0000000 
-0.6727198 
0.1878524 
-10.0146337 


HIGHEST 

SIDELOBE 

-13 dB 

-43 dB 

-72 dB 

-98 dB 


3-dB WIDTH 

0*89 DF 

1.30 DF 

1.61 DF 

1.86 DF 


6-dB WIDTH 

1.20 DF 

1.81 DF 

2.25 DF 

2.62 DF 


! for OUR use I THINK THE 3RD-ORDER WINDOW IS THE MOST APPROPRIATE, AND 
JtHE PROGRAM USES THIS WINDOW FUNCTION. SINCE W3 - 0 IN THIS CASE, THE 
ACTUAL CONVOLUTION CONTAINS ONLY 3 TERMS. SINCE W0 - 1 , THE LAST TERM 
iIS OUST X<j), WITHOUT A MULTIPLICATION. THUS THE WINDOWING OF EACH 
iPOINT INCLUDES 5 ADDITIONS AND 2 MULTIPLICATIONS FOR EACH QE REAL 

jAND IMAGINARY COMPONENTS. IF ONLY VERY SHORT BLOCKS (SAY NX ■ 64) ARE 
tTO BE PROCESSED, THE 2ND-ORDER WINDOW MIGHT BE MORE APPROPRIATE SINCE 
tTHE FILTERS ARE A BIT NARROWER IN THAT CASE. 
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5 V PROGRAM 

! * 

,****************»******** ******* ************************* 
5 

! DEFINITIONS : 


* 
* 

*************** 


1000 

DEFINT I-N 

1010 

NX = 1024 

1020 

NXM1 = NX-1 

1030 

DIM X( NX > 

1040 

DT - 0.01 

1050 

DF - 1 . /( DT*NX ) 

1060 

NMAX - 30. /DF 

1100 

NF - NX\2 

1110 

MF - LOG< NF )/LO 

1120 

NFD2 *= NF\2 

1130 

NFM1 - NF-1 

1140 

DIM CA(NF) 

1150 

DIM SA(NF) 

1160 

PI2 « 6.2831853 

1200 

W1 - -0.5859613 

1210 

W2 » 0.0922278 


'VARIABLES STARTING "l"-«N" ARE INTEGERS 
’DATA BLOCK SIZE 


’INPUT DATA ARRAY X( 1 >, 1 = 0 

'INPUT SAMPLE TIME IN SECONDS 
'SPECTRAL LINE SPACING IN Hz 
•NUMBER OF HIGHEST SPECTRAL LINE <- 


NX-1 


30 Hz 


'FFT SIZE (NR OF COMPLEX POINTS) 
'NF - 2~MF 


'COSINE AND SINE ARRAYS, 

• ANGLES FROM 0 TO 1 80*( NX- 1 )/NX 


'WINDOW COEFFICIENTS 


,*****************************„ ****************************************** 

5 

; CALCULATE SINE AND COSINE ARRAYS: 

: 

! 5 0NE 0NLV 0NCE * NOTE THAT THE FFT ROUTINE USED SINES AND 
{COSINES FROM 0 TO <180 DEGREES IN STEPS OF 360/NF DEGREES AND THE RFAf 
{TRANSFORM EXTRACTION ROUTINE USES SINES AND COSINES UP To ’<90 DEGREES^ 
J60/NX DEGREES * THUS THE ODD NUMBERED POINTS IN T^SE SrSs 
ANGLES ABOVE 90 DEGREES ARE NOT USED. IF MEMORY SPACE IS A PROBLEM 
BE CONDENSED AND ONLY POINTS FROM 0 TO <180 DEGREES IN * 
;STEPS OF 360/NF CALCULATED FOR USE BY THE FFT, IN THIS CASE THF v/ai iifc 

; NFFnpn V T ?t 5c A a, TRANSF0RM EXTRACTION WILL HAVE TO BE 1 CALCULATE D^^^^ 
;NEEDED. IT IS ALSO POSSIBLE TO USE A TABLE CONTAINING VALUES FOR anpi f c 

? Mnpc r^Di frISc 45 DEGREES « BUT THIS MAKES THE LOOKUP PROCEDURE SOMEWHAT 
;MORE COMPLICATED, INVOLVING A POSSIBLE CHANGE OF SIGN OR EXCHANGE OF qtm 
{AND COS VALUES DEPENDING ON THE SEMIQUADRANT OF THE ANGLe! SIN 


1500 DA - PI 2/NX 

1510 FOR 0 » 0 TO NFM1 

1520 A = DA*J 

1530 CA( J ) = COS(A) 

1540 SA( 0 ) *> SIN( A > 

1550 NEXT J 


'ANGLE INCREMENT 
•CALCULATE SINES AND COSINES 



» 


{SAMPLE INPUT DATA INTO ARRAY X: 


{THIS CODE MUST BE WRITTEN TO SAMPLE THE INPUT SIGNAL EVERY 0 01 SFrnnin 
{AND STUFF SUCCESSIVE VALUES INTO X. ONLY NX VALUES NEED BE TA^FN 
{ENTER HERE TO START PROCESSING A NEW BLOCK OF DATA. 


{2000 CODE YET TO BE WRITTEN 
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* B IT REVERSALS 

* THE REAL COMPACT SAMPLES IN X ARE NOW TREATED AS A COMPLEX ARRAY OF NF 
(REAL, IMAGINARY) PAIRS. PAIRS AT BIT-REVERSED ADDRESSES ARE NOW SWAPPED 
PRIOR TO THE FFT ROUTINE. 


; THE FOLLOWING ROUTINE USES LOGICAL OPERATORS. 


3000 

3010 

3020 

3030 

3040 

3050 

3060 

3070 

3080 

3090 

3100 

3110 

3120 

3130 

3140 

3150 

3180 

3190 


K - 1 'INITIAL BIT-REVERSED ADDRESS 

FOR 0 - 1 TO NFM1 'SCAN ALL ADDRESSES EXCEPT THE FIRST 

L - NFD2 'INITIAL BIT MASK 

IF NOT( K AND L) GOTO 3070 'THIS BIT SET? 


K - K AND NOT L 
L - L\2 
GOTO 3030 
K * K OR L 
IF K <- J GOTO 3190 
02 - U+J 
K2 - K+K 
T - X( 02 ) 

X( 02 ) - X(K2> 

X(K2) - T 
T - X< 02+1 ) 

XI02+1 ) - X( K2+1 > 

X( K2+1 > - T 
NEXT 0 


• YES - CLEAR IT, 

' SHIFT MASK RIGHT, 

' AND CONTINUE 

•BIT WAS CLEAR, SET IT 

■SKIP IF ADDRESS SAME OR ALREADY DONE 

•ARRAY ADDRESS OF LOWER REAL ^COMPONENT 

;ARRAY ADDRESS OF UPPER REAL COMPONENT 

■SWAP REALS 


'SWAP IMAGINARIES 


j THE FOLLOWING VERSION USES ONLY INTEGER ARITHMETIC BUT IS A BIT LESS 
; EFFICIENT . 


; 3000 

FOR 0 - 1 TO NFM1 

; 3010 

K1 - 0 

; 3020 

K - 0 

i 3030 

FOR I - 1 TO MF 

; 3J0T4J0T 

K2 - K 1 \2 

*3050 

K » ( K-K2 ) *2 + 1 

; 3060 

K 1 - K2 

; 3J0T7J0T 

NEXT I 

;308 0 

IF K <« J GOTO 

; 3090 

02 * 0+0 

i 3100 

K2 - K+K 

; 3 1 10 

T - X ( 0 2 > 

; 3120 

X< 02 ) - X(K2> 

; 3 1 30 

X ( K2 ) - T 

; 3 1 40 

T - X < 0 2+1 > 

; 3 1 50 

X ( 02+ 1 ) - X( K2 + 

; 3 1 80 

X ( K2+ 1 ) - T 

; 3 1 90 

NEXT 0 


'SCAN ALL ADDRESSES EXCEPT THE FIRST 
'INITIALIZE DIVISOR 
'INITIALIZE BIT-REVERSED ADDRESS 
sBIT-REVERSE ALL MF BITS 


•SKIP IF ADDRESS SAME OR ALREADY DONE 
■ARRAY ADDRESS OF LOWER REAL COMPONENT 
jARRAY ADDRESS OF UPPER REAL COMPONENT 
•SWAP REALS 


'SWAP IMAGINARIES 
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1/D8/J 


FFT ROUTINE: 

THIS ROUTINE USED DECIMATION IN TIME. THE COMPLEX DATA INPUT IS IN BIT- 
REVERSED ORDER. 

1. THE OUTER LOOP MAKES MF PASSES THROUGH THE DATA. THE INITIAL 
NUMBER OF SUB-DFT'S IS 1, AND DOUBLES ON EACH PASS. THE INITIAL 
INCREMENT IN ANGLE BETWEEN SUB-DFT'S IS 180 DEGREES, ANO HALVES 
ON EACH PASS. 

2. THE MIDDLE LOOP PROCESSES THE SUB-DFT'S. AT EACH SUB-DFT 
ANGLES INCREASE IN NATURAL ORDER AS 0 <1-1 >* 0, 90 <I-2>; 0, 45, 
90, 135 < I -3 > ; 0, 22.5, 45, ... <I»4); AND SO ON. 

3. THE INNER LOOP CALCULATES THE BUTTERFLIES FOR A GIVEN ANGLE 
IN EACH SUB-DFT. EACH DEC I MATION- I N-T I ME BUTTERFLY HAS THE FORM 

Re< j ) <- Re( j ) + E Re( k )*C+Im(k)*S] 

I m< j ) <- Im(j ) + EIm(k)*C-Re<k>*S] 

Re< k ) <- Re< j )-CRe(k)*C+Im<k)*S] 

Im(k) <- I m< j )- C Im< k ) *C-Re< k )*S] 

WHERE C AND S ARE THE COSINE AND SINE OF THE ANGLE. 


3200 

LM » 1 

•INITIAL NUMBER OF SUB-DFT'S - 1 

3210 

3220 

LM2 « LM+LM 
IDA - NF 

•INITIAL ANGLE INCREMENT 

» 180 DEG 

3230 

FOR I - 1 TO MF 

'OUTER LOOP 


3240 

IA - 0 

•INITIAL SUB-DFT ANGLE * 

0 DEGREES 

3250 

FOR L - 1 TO LM 

'MIDDLE LOOP 


3260 

C - CA< IA) 

'LOOK UP COS AND SIN FOR 

THIS ANGLE 

3270 

3280 

S - SA< IA) 

FOR U-L-l TO NFM1 STEP LM2 ; I NNER LOOP 


3290 

OR - 0+0 

•ARRAY ADDRESS OF LOWER 

REAL COMPONENT 

3300 

01 - OR+1 

‘ARRAY ADDRESS OF LOWER 

IMAG COMPONENT 

3310 

KR » 0R+LM2 

•ARRAY ADDRESS OF UPPER 

REAL COMPONENT 

3320 

KI - KR+1 

•ARRAY ADDRESS OF UPPER 

IMAG COMPONENT 

3330 

TR - X ( KR > *C + X < K I ) 

*S 'BUTTERFLY 


3340 

T I - X < K I >*C - X< KR ) 

*S 


3350 

3360 

3370 

3380 

3390 

X ( KR ) - X< OR )-TR 
X ( K I ) - X<01 )-T I 
X < 0 R > « X< OR )+TR 
X < 0 I ) » XCOI > +T I 
NEXT 0 

'END OF INNER LOOP 


3 400 

IA - IA+IDA 

'BUMP ANGLE 


3410 

NEXT L 

'END OF MIDDLE LOOP 


3340 

LM = LM2 

•DOUBLE NUMBER OF SUB-DFT’S 

3430 

3440 

LM2 - LM2+LM2 
IDA - I DA\2 

•HALVE ANGLE INCREMENT 


3450 

NEXT I 

•BOTTOM OF OUTER LOOP 
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J 

; EXTRACT CLOSE-PACKED REAL DATA TRANSFORM: 

!now WE MUST EXTRACT THE SPECTRUM OF THE INPUT DATA AS IF IT HAD BEEN 
\ CALCULATED FROM A SEQUENCE OF NX COMPLEX INPUT VALUES ALL OF WHOSE 
t IMAGINARY COMPENENTS WERE ZERO. THE ALGORITHM FOR THIS IS AS FOLLOWS: 


; 


i 

1. FOR POINTS j - 1 

, 2, . NF/2-1 CALCULATE 

? 

? 

R©( j ) <- R©< j > + Re<k)-tRe( j >-R©(k)3*S+CIm< j > + Im(k)3*C 


Im( j ) <- I m< j )- Im< k >- CR©( j >-Re< k >3 *C- C Im< j > + Im< k >3 *S 

; 

Re( k > <- Re( j >+Re< k >+ CR©< j >-Re< k > 3 *S- C Im( j > + Im<k)3*C 

; 

Im(k) <- - Im( j >+Im<k)-CRe(j >-R©< k ) 3 *C- C I m< j ) + Im( k ) 3 *S 

. 

WHERE k - NF -j, A 

- j *360/NX DEGREES , C » COS(A) l AND 

; 

S - SI N( A > . 


; 

? 

WE ALSO HAVE THE SPECIAL POINTS: 

i 

; 

R©( 0 > <- CRe< 0> + Im< 0 >3 

X 

Im( 0 ) <- 0 


; 

i 

Re( NF/2 ) <- 

R©( NF/2 > *2 

? 

Im< NF/2 ) <- 

- Im( NF/2 >*2 

» 

! 

AND (THOUGH THESE VALUES ARE NOT CALCULATED): 

i 

; 

Re(NF) <- CR©( 0 >- Im< 0 ) 3 *2 

: 

Im< NF ) IS UNKNOWN 

3500 

X(0> - X( 0 )*2 . 

'SPECIAL CASE FOR X(0> 

3510 

X( 1 ) - 0. 


3520 

FOR U - 1 TO NFD2-1 

'LOOP FOR GENERAL CASE 

3530 

UR - U+U 

•ARRAY ADDRESS OF LOWER REAL COMPONENT 

3540 

UI » UR+1 

'ARRAY ADDRESS OF LOWER 1MAG COMPONENT 

3550 

KR - NX-UR 

'ARRAY ADDRESS OF UPPER REAL COMPONENT 

3560 

K I - KR+1 

•ARRAY ADDRESS OF UPPER IMAG COMPONENT 

3570 

C - CA( 0 ) 

'LOOK UP COS AND SIN FOR THIS ANGLE 

3580 

S - SA( J ) 


3590 

RS - X( UR )+X( KR ) 

'REAL SUM 

3600 

RD - X ( UR ) -X ( KR ) 

'REAL DIFF 

3610 

QS « X(UI ) +X ( K I > 

' IMAG SUM 

3620 

QD « X(UI ) — X { K I > 

‘ IMAG DIFF 

3630 

TR * QS*C-RD*S 

■TERM FOR REALS 

3640 

TI - RD*C+QS*S 

•TERM FOR IMAGS 

3650 

X( UR > - RS + TR 


3660 

X(UI > - QS-TI 


3670 

X( KR > - RS-TR 


3680 

X ( KI > - -QS-TI 


3690 

NEXT U 

■BOTTOM OF LOOP 

3700 

X ( NF > - X ( NF ) *2 ♦ 

■SPECIAL CASE FOR X(NF/2) 

3710 

X ( NF + 1 > - X ( NF + 1 )*2 . 
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;WIND0W SPECTRAL DATA TO IMPROVE FILTER SHAPE: 

• W I NDOW^F UNCTION ! 0F WIND0WING AB0VE * ™ IS ROUTINE IS FOR A 3RD-ORDER 


3800 

3810 

3820 

3830 

3840 

3850 

3860 

3870 

3880 

3890 

3900 

3910 

3920 

3930 

3940 

3950 

3960 


'INITIALIZE PREVIOUS VALUE ARRAY , Re< 2 ) 
Re( 1 > 


R2 - X( 4 > 

R1 ■ X( 2 ) 

R0 = X(0> 

Q2 =■ -X( 5 > 

Q1 = -X< 3 > 

Q0 ■ — X < 1 > 

F °ii J0 T ° NMAX * 2 STEP 2 'WINDOW POINTS <« 30 Hz 

» >-wiic SK5f aS:*5WSiSa P ?5Kr component 

X(un - <Q2 + X(O1+4))*W2 + (Q1+X(O1+2))*W14.Q0 • IMAG COMPONENT 


Re( 0 ) 

I m( 2 ) 

I m{ 1 ) 

I m< 0 ) - 


0 . 


R2 
R 1 
R0 
Q2 
Q1 
00 
NEXT 


R1 
R0 

X ( 00+2 ) 

Q1 

Q0 

X< Jl+2) 


'SHIFT PREVIOUS ARRAY 
AND INSERT NEW VALUE 


'END OF LOOP 


******************* 


********************* 




CALCULATE MAGNITUDES: 

THE MAGNITUDE OF EACH SPECTRAL VECTOR IS CALCULATED AS 

MAGNITUDE! j ) * SQRTCRe! j )*Re(j > + Im(j )* Im< j ) J *SF 

WHERE SF - 1/NX IS A SCALING FACTOR, USED FOR 2 REASONS: 

i,*,- SCALING FACTOR SETS THE PROCESSING GAIN WITH <:f * i /mw 

THE GAIN FOR A SINGLE-COMPONENT SIGNAL IS 1. THAT is" IF tuf^*’ 
INPUT SIGNAL IS A SINE-WAVE AT FREQUENCY n*DF WITH PEAK AMPLlTiinr 
1.0, THE SPECTRUM AT FREQUENCY n*DF WILL ALSO BE 1.0 AMPLITUDE 

THE BL0CK C ^l7? G N v AC In R Tur K cLI HE PROCESSING GAIN INDEPENDENT OF 
nrrrrn?^ SIZE NX, SO THE SAME SIGNAL CAN BE PROCESSED WTTH 

?J E Lv e ^s; E ?„ E F °s,^^ri/?;r f s ,^T 3 w;- 


4000 

4010 

4020 

4030 

4040 

4050 


cno " 1 • 'i* ‘ SCAL E FACTOR 

F0 ^° “ J 6 T0 NMAX 'PROCESS POINTS O 30 Hz 

■ *** 'ADDRESS OF REAL COMPONENT 

w7nT 0Z t l 'ADDRESS OF IMAG COMPONENT 

neJU ■ ^T<»"*>**<«>gmj*x<K2>£SF ' calculate magnitude 


: END OF PROGRAM. 
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TU r dcct \jay TO TEST THE PROGRAM IS TO GENERATE A SYNTHETIC INPUT ARR ^Y 
OF DATA POINTS CONTAINING SIGNALS AT KNOWN FREQUENCIES, AND; PROCESS THEM 
AND CHECK THAT THE RESULTING SPECTRUM IS CORRECT. TO GENERATE A SIGNAL 
ARRAY WHICH CONTAINS A SIGNAL AT A FREQUENCY Q*DF USE THE FALLOWING CODE 
BEFORE THE BIT-REVERSAL ROUTINE: 


2000 Q - . . . 

2010 DP - P I 2*Q/NX 

2020 FOR I - 0 TO NXM1 
2030 X< I > ■ SIN<DP*I ) 

2040 NEXT I 

t Q IS A REAL NUMBER, AND SHOULD 
5 THAT THE ARRAY X WILL CONTAIN 


'FREQUENCY IN UNITS OF DF 
•PHASE INCREMENT PER SAMPLE 
'DATA GENERATION LOOP 
■CALCULATE A DATA POINT 
'END OF LOOP 

BE CHOSEN SO THAT 0 <- Q <* NX/2. 
I CYCLES OF A SINEWAVE SIGNAL. 


NOTE 


actfr PROCESSING THE SPECTRUM SHOULD CONTAIN ONLY COMPONENTS NEAR THE 
FREQUENCY Q*DF. THAT IS. SPECTRAL LINES XC j > SHOULD BE SMALL EXCEPT FOR 
NEAR Q if Q IS AN INTEGER, THE SPECTRAL LINE X(Q) SHOULD BE 1 AND 
ALL THE OTHERS SHOULD BE 0 (OR OF ORDER IE-6, DEPENDING ON THE ROUND-OFF 
FRRORS IN THE PROCESSING)* USING NON-INTEGER VALUES OF Q WILL RESULT IN 
A SPECTRUM WITH A CLUSTER OF RESPONSES NEAR X(INT(Q>> BUT L ; I N E S FAR AWAY 

frSm E 0 Should be small the falloff in response in the spectrum as one 
JiSSes aS2? f£om the lime X(INT<Q>> depends on the effectiveness of the 

WINDOW FUNCTION, AND THE SHAPE OF THE SYNTHESIZED FILTER RESPONSE CAN BE 
PLOTTED BY USING VARIOUS VALUES OF Q. 

NOTE THAT IF NMAX < Q <* NX/2, THE SIGNAL GENERATED LIES I H UPPER 

PART OF THE SPECTRUM, WHICH IS THROWN AWAY (DON'T EXAMINE SPECTRAL LINES 
y/ . \ uhfrF i > NMAX AS THEY HAVEN'T BEEN CALCULATED)* ALSjO, IF 
NX/2 < Q < NX, THE SIGNAL WILL BE ALIASED IN THE SPECTRUM, APPEARING AT 
A FREQUENCY OF NX/2-Q. 

*********************************************************** ************* 


j 
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ULF TREE POTENTIALS AND GEOMAGNETIC PULSATIONS (U) 
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(the numbers are 

9.XRAY >= M5 





transmitted with 

10.XRAY >= 

XI 





the alerts) 

11. A INDEX 

>= 

20 

(B0U) 

(forecast or observed) 


12. A INDEX 

>= 

30 

(B0U) 

fl II If 



13. A INDEX 

>= 

50 

(B0U) 

M If It 



1 4.K INDEX 

>= 

4 

(B0U) 

(fcst, obsvd,or in 

progress) 


15.K INDEX 

>= 

5 

(B0U) 

h it ii ii 

ii 


1 6.K INDEX 

>= 

6 

(B0U) 

u if ii h 

ii 


17.K INDEX 

>= 

5-5 

(B0U) 

if ii if n 

if 


18. SUDDEN STORM COMMENCEMENT (forecast or observed) 

1 9.STRATWARM 

20.SST ALERT: RADIATION >= 10 MILL I REMS (forecast or 
observed) 

LINE SEVEN... MESSAGES 


where I 1 1 1 !••«* unformatted message with up to 60 characters wl I I 

sometimes be broadcast In this line. Messages that 
will take more than one line will be broadcast 
sequentially over several minutes. These sequential 
messages will start with MMMM at the beginning of 
the first line. NNNN at the beginning of a line 
will Indicate that the message is complete. 

In the future routine products will be transmitted 


using this line at preselected times, 
schedule of these products follows: 

A tentative 

PRODUCT 

TIME 

Current message schedule 

OOZ 

SMM (Solar Maximum Mission) 
Observing Plan 

01Z,07Z,13Z,19Z 

Region Report (body of 0030Z report) 

02Z,08Z, 1 4Z,20Z 

27 Day Forecast 

03Z,09Z,15Z,21Z 

SGAS 


(Solar Geophysical Activity Summary) 
(body of 0245Z report) 

04Z,10Z,16Z,22Z 

SDF (Report of Solar and Geophysical 
Activity. .former ly Space Disturbance 
Forecast) (text portion of 2200Z 
report) 

05Z,11Z,17Z,23Z 

GEOALERT (body of 0330Z GEOALERT) 

06Z,12Z,18Z 
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LINE FIVE FORECASTS AND WARNINGS (continued) 


Minute 4 

I I 1 1 1 1 I II 1 1 I ... : 

SOLQU I ET 

where 11111...= Description of forecast solar activity 

This will usually be: which Indicates 

SOLQU I ET only sporadic weak solar activity, 

(class C or smaller flares) 

SOLALERT JJ/KK M Class flarqs, or an 

i ncrease I n NJ I eve I 
activity fronji days JJ 
to KK 

MAJOR FLARE ALERT JJ/KK QXXYY 

region located at QXXYY 
(see minute 3) has a 1 5% 
or greater prob of X class 
flares from days JJ to KK 

SOLNIL end of active period or begin Ing 

of period of ! very low activity 
(used to end alerts) 


LINE SIX ALERTS 

(rotates through all activated alerts) 

1 1 1 1 1 II 1 1 . : 

ALERT NIL 

where 11111...= descriptive text of alert conditions that are met 
These are The current alerts 
0. ALERT NIL (when no alerts are turned on) 

1.245 MHZ BURST > 100 SFU 

2.245 MHZ RADIO NOISE STORM > 5 TIMES BKGND 
3.10 CM FLARE > 100$ ABOVE BKGND 

4. TYPE II RADIO BURST 

5. TYPE IV RADIO BURST 

6. SUSPECTED PROTON FLARE 

7. PROTON EVENT 10 PART/CM2*SEC*STER > 10 MEV 
8 PROTON EVENT 100 PART/CM2*SEC*STER > 100 MEV 
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MAGALERT JJ/KK or Increase In or continued 

MAGALERT MINOR JJ/KK high magnetic activity 

(JJ/KK are the first and 
last calendar days of 
activity) 


MAGNIL end of active period or beginning 

of period of very low activity 
(used to end alerts) 


Minute 3 

1 1 I I I 1 1 II I I I ... 5 

PROTONN I L 

where Mill...- Description of forecast proton activity 

This will usually be: which Indicates: 

PROTONQUIET no proton enhancements measured 

in, or forecast for the near 
earth environment. 

PROTON FLARE ALERT JJ/KK QXXYY 

Protons expected near earth 
because the solar region 
located in quadrant Q (1=NE, 

2=SE, 3-SW, 4=NW) , Central 
meridian distance XX, latitude 
YY, has a 15? or greater 
probabilty of a proton 
flare from days JJ to KK 

PROTON ARRIVAL ALERT PP/FF HHmm 

Protons expected near earth 
on day PP from flare which 
occurred on day FF at HHmm 

PROTONALERT JJ/KK Protons expected to meet the 

SESC alert criteria (see line 6 
ALERTS) in the near earth 
environment from days JJ to KK 

PROTONN I L end of active period or begtning 

of period of very low activity 
(used to end alerts) 
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LINE FIVE.... FORE CASTS AND WARNINGS 
Minute 0 

doy-doy Mnn/nn/nn Xnn/nn/nn Pnn/nn/nn lOCMrrr/rrr/rrr/rrr 

348-350 M01/01/01 X01/01/01 P01/01/01 10CM1 02/100/098/*** 

where doy-doy = day of year, first and last days of the forecast 
nn/nn/nn = percent probability for occurrence of activity 

on the first, second, and last day of the forecast. 
These apply to M flares, X flares, or Proton events 
as indicated. Example 01/01/01 

rrr/rrr/rrr/rrr = forecast 10 CM solar radio flux for the 3 forecast 

days, such as 102/100/098/*** The last group of 
rrr is for the last doy plus one and will be *** 
until approx 18z on the first doy. At this time 
values will be added. 

Minute 1 

doy-doy AFaaa/aaa/aaa APaaa/aaa/aaa Kkkkk kkkk kkkk kkkk : 

348-350 AF01 8/ 015/015 AP020/020/020 K**** **33 3444 3333 

day of year, first and last days of the forecast 
forecast A values for the forecast per iod, either 
AF, or AP as indicated. Example 018/015/015, or 
20/20/20 

Forecast KP values (without the 0,+,or-). 

The first of these 16 consecutive K values 
is for 00Z through 03Z on the first Idoy listed. 

Once a forecast K value has been measured it will 
appear in line two and the forecast va I ue in this 
line will be replaced by an *. Example ***# ##*3 
3444 3334. At approximately 18Z the remaining 
values will be updated with a new forecast. 


Minute 2 

I I I 1 1 1 1 1 1 1 1 ... : 

MAGALERT MINOR 

where I INI...® Description of forecast magnetic activity 
This will usually be: which indicates: 

MAGQUIET only sporadic weak 

geomagnetic activity 
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Minute 3 


doy 10CM rrr 90 DAY MEAN mmm SSN sss XRAY BGD lx.x 
346 10CM 104 90 DAY MEAN 110 SSN 090 XRAY BGD B2.5 


where 


doy 

rrr 

mmm 

sss 

lx.x 


= day of year (yesterday). Such as 346 
= 10 CM solar radio flux In Solar Flux Units, 
for 1700Z on doy Indicated, such as 104 
= 90 day mean of the 10 CM Solar radio flux 
through yesterday, such as 110 
= Sun Spot Number for yesterday, such as 090 
- Xray Background level for yesterday expressed on 
the B,C,M,X scale, where B = 10“7 watts / meter2 

£ B 10” , M = 10“*, and X = 10”^. In the example, 
B2.5 = 2.5 * 10-7 wat+s / me ter2. 


Minute 4 


doy G0ES5 (W075) FLUENCE 1MEV f.fE+f 10 MEV f.fEtf 
346 G0ES5 (W075) FLUENCE 1MEV 2.8E+4 10 MEV 5.2E+3 


where doy = day of year (yesterday), such as 346 

f.fE+f * dally total fluence (protons / (cm2*clay*ster ) ) 

greater than 1MEV or 10MEV (as indicated) expressed 
as an exponential number, such as 2.8E+4, or 5.2E+3 


LINE FOUR.... HOURLY IONOSPHERIC DATA 

Each 24 hour set of data will take 3 minutes to transmit 
(8 hours of data per line). EMPTY Indicates that no data 
is being transmitted. 

sss ddddd doy hh nn.n nn.n nn.n nn.n nn.n nn.n nn.n nn.n : 

BOU TEC 347 16 09.5 14.0 20.3 22.8 22.3 22.8 **** **** 


where sss 
ddddd 


hh 

nn.n 


* station identifier such as BOU: 

= data Identifier, such as TEC (Total Electron 
Count), F0F2 (Critical Frequency of the F2 Layer), 
M3000 (Optimum Frequency for a Single Hop 3000 KM 
Transmission), FMIN (Minimum Usable Frequency), 
or FOES (Critical Frequency of the Sporadic 
E Layer). 

- hour of first data group, such as 16 
= data group. There are eight consecutive data 
groups. Missing data, or data not yet 
determined will have *#**. Data for TEC 
is In units of 10+16 electrons / meter2, 
for the others It is MHZ. Example 09.5 


42 

Approved For ReleasebtMGi(ASSJIHGD6 -00788R001 8001 90001 -3 


Approved For Release 2(|^^^C^^|^^)0788R001 8001 90001 -3 


Minute 1 


doy hh-hh K k/ggg BO Knnnn nnnn PL Knnnn nnnn 10CM rrr 
347 18-21 K 3/033 BO K3334 333* PL K4334 433* 10CM 104 


where doy 
hh-hh 

k/ggg 

nnnn nnnn 


rrr 


= day of year (current day), such as 347 
= beginning and ending hours (in Universal 
Time or Z) for the 3 hr period of tijie K 
index to follow, such as 18-21 
= Boulder K value / Gammas of deflection from the 
quiet day curve, such as 3/033 
= K indices determined so far for this day (Boulder 
or planetary). The first value is for 00Z-03Z, 
and each index is for the 3 Hr period following 
the K index preceding it. Missing values will 
have an *. Example 3324 333* or 4334 433* 

= 10 CM solar radio flux in Solar Flux Units, 

(10”22 watts / (meter2*Hz)) as measured by 
Ottawa at approximately 1700Z each day. This value i 
for the doy Indicated at the start of the line. 

Until this value Is available *** will be broadcast. 
Example 104 


Minute 2 

doy B0 Aaaa Knnnn nnnn PL Aaaa Knnnn nnnn TEDttt : 

346 B0 A0 17 K3334 4323 PL A024 K4333 4434 TED*** 

where doy = day of year (yesterday), such as 346 

aaa = A index (Boulder or Planetary) for yesterday, 
such as 017, or 024 

nnnn nnnn = K indices (Boulder or Planetary) for yesterday, 

(starting at 00-03Z) missing values will contain *. 
Example 3334 4323 or 4333 4434 
ttt = Total Energy Deposition, an estimate of the total 
power being dissipated in one hemisphere's auroral 
atmosphere by precipitating particles. The units 
are gigawatts. This is a daily average for the doy 
Indicated. Estimates are made using particle data 
from Tiros satellites. Daily averages vary from 
near 005 on a quiet day up to about 100 on a very 
active day. If this data is unavailable *** will 
be used. Example *** 
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LIME TWO GOES GEOSYNCHRONOUS MAGNETIC DATA 

(1 minute values broadcast each minute) 

G0ES5 HPsgggg HEsgggg HNsgggg G0ES6 HPsgggg HEsgggg HNsgggg : 

G0ES5 HP+0123 HE-0023 HN+0033 G0ES6 HP+0123 HE-0023 HN+0033 

where sgggg = sign (+ for positive, - for negative) and value 
(In gammas) of the magnetic field components at 
geosynchronous altitudes (approximately 22,000 
miles). These are measured by the GOES 5 (W075) 
and GOES 6 (W135) satellites. These components are 
HP (parallel to the Earth's rotation ax!s,+ is 
north), HE (directed earthward, + is down) and HN 
(perpendicular to the other two, le east or west, 

+ is west). Example +0123. 


LINE THREE.... INDICES 
Minute 0 

1MEV p.pE+p 10MEV p.pE+p 100MEV p.pE+p PCA aa.aDB NTRN nnnn : 

1MEV 3.1E+0 10MEV 2.5E+0 100MEV 1.2E+0 PCA 00.0DB NTRN 0451 

where p.pEp = number of protons / (sec*cm2* s ter) 

at greater than 1MEV, 10MEV,or 100MEV energies 
(5 min avg), such as 3.1E+0 
aa.a = DB of polar cap absorption at 30 MHZ as 

measured at Thule, (15 Min value), such as 00.0 
nnnn = Thule neutron monitor data In millivolts (15 min 
value). The background level is 460 to 470 mv. 
Significant deviations are + or - 10$ or less than 
410 mv and greater than 510 mv. A sharp increase 
that lasts for a few tens of minutes to hours is a 
Ground Level Event that Is caused by relativistic 
particles associated with solar flares. A sig- 
nificant decrease that lasts for several hours 
to several days Is a Forbush Decrease. This 
decrease is the scattering of background cosmic 
rays due to plasma density enhancements and 
associated magnetic field enhancements in the near 
earth solar wind. These enhancements may or may 
not be solar flare related. 
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Appendix C 

FORMAT FOR SATELLITE BROADCAST OF SPACE ENVIRONMENT SERVICES (U) 

(15 February 1984) 


There are seven sixty character lines. Each is followed by a carriage 
return and a line feed. 

1 234567890 1 234567890 1 234567890 1 234567890 1 234567890 1 234567890 

In the examples of lines, CAPITAL letters are caption Information 
that appears as it does in the broadcast. Small letters represent 
data that are explained below. The : indicates the 61st! space 
(: Is not broadcast). Any data not available for broadcast will be 
replaced by **#. 


LINE ONE - GOES X-RAY AND TOTAL MAGNETIC FIELD 
(1 min values broadcast each minute) 

! 

doy hhmm G0ES6 lx.x G0ES5 lx.x BOU TF ggggg GAMMAS : 

347 2153 G0ES6 B2.0 G0ES5 B2.0 BOU TF 55397 GAMMAS 

I 

where doy = day of year, such as 347 

hhmm = hours and minutes, such as 2153 

lx.x = x-ray level, (1 min avg), expressed with the B,C,M,X 
scale, where B = 10-7 wa++s / me +er2, C = 1 0“6 
M = 10~5 anc | x = 10 -4. 

X-rays are measured in a 1-8 angstrom channel on the 
GOES geosynchronous satellites. An increase in X-ray 
levels Indicates solar flare activity. 

In the example B2.0 = 2.0 * 10"7 watts / meter2. 

ggggg = Boulder mag field in gammas, (1 min avg), 
such as 55397 
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FORMAT FOR SATELLITE BROADCAST OF SPACE ENVIRONMENT SERVICES (U) 
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desirable, because the tree ring data could provide a unique 
record of past ULF and higher-frequency geomagnetic 
activity. 
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